New Phases of Water Ice Predicted at Megabar Pressures 
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Based on density functional calculations we predict water ice to attain two new crystal structures with Pbca 
and Cmcm symmetry at 7.6 and 15.5 Mbar, respectively. The known high pressure ice phases VII, VIII, X, 
and Pbcm as well as the Pbca phase are all insulating and composed of two interpenetrating hydrogen bonded 
networks, but the Cmcm structure is metallic and consists of corrugated sheets of H and atoms. The H atoms 
are squeezed into octahedral positions between next-nearest O atoms while they occupy tetrahedral positions 
between nearest atoms in the ice X, Pbcm, and Pbca phases. 

PACS numbers: 



Water ice is one of the most prevalent substances in the so- 
lar system, with the majority of it existing at high pressures in 
the interiors of giant planets OJJ] . Uranus and Neptune are as- 
sumed to consist largely of a mixture of water, ammonia and 
methane ices at pressures up to 8 Mbar [2], while Jupiter and 
Saturn almost certainly have large dense cores consisting of 
differentiated rock and ice components at pressures on the or- 
der of 10 Mbar in Saturn |Q] and 39-64 Mbar in Jupiter |H]. 
The behavior of water ice under these extreme conditions is 
not yet well understood because static diamond anvil cell ex- 
periments have not yet reached beyond 2.1 megabars 0-01 • 
Shock wave experiments [8] have reached higher pressures 
but they also heat the sample significantly so that it melts 
for the highest pressures. However, dynamic ramp compres- 
sion techniques [9, 10] are expected to reach high pressures at 
comparatively low temperatures, where so far only theoretical 
methods have predicted the state of water ice. 

The known phase diagram of water is extremely rich, with 
at least fifteen forms of solid ice observed experimentally ill ill 
and one high-pressure phase predicted theoretically 

El. The 

high-pressure region of the phase diagram is comparatively 
simple. The molecular ice VIII phase, which forms at low 
temperatures for pressures above ~15 kbar, consists of a bcc 
array of oxygen atoms with an ordered arrangement of hydro- 
gen atoms arranged along the tetrahedral directions, bonded 
with a short covalent bond to one oxygen and a longer hydro- 
gen bond to the other. Increasing pressure to 0.7 Mbar results 
in the symmetrization of these bonds 11311 . with the distinction 
between covalent and hydrogen bonds being lost as the hy- 
drogen occupies the midpoint between the two oxygens. This 
is referred to as phase X, the highest-pressure phase that has 
been observed experimentally. In 1996 Benoit et al. lfl2ll used 
density functional theory (DFT) to predict a higher-pressure 
phase of ice with Pbcm symmetry and twelve atoms per unit 
cell to become stable at approximately 3 Mbar. It was recently 



shown within DFT that this pha se forms by a dynamic insta- 
bility in the ice X lattice ljl4ll5ll . 

At temperatures in excess of approximately 2000 K, high- 
pressure ice transitions to a superionic phase 1 3- 18 ] in which 
the hydrogen atoms become mobile while the oxygen atoms 
do not, while at higher temperatures still the oxygen atoms 
also become mobile and the entire structure melts 1 191 In the 



recent work by French et al. 112011 it was shown that the bcc 
structure of the oxygen lattice in the superionic phase appears 
to be maintained up to very high pressures. This work also in- 
dicated that at densities below approximately 5 g cm~ 3 , corre- 
sponding to pressures around 6 Mbar, the hydrogen atoms are 
found to be strongly associated with the four tetrahedral sites 
surrounding each oxygen as in ice X, but at higher pressures 
the hydrogen atoms show an increasing preference for the six 
octahedral sites surrounding around each oxygen. It is thus 
natural to ask whether the occupation of these sites may lead 
to the formation of a novel crystalline phase of ice that would 
be more stable than the Pbcm phase at high pressure. 

To investigate the plausibility of the occupation of the oc- 
tahedral sites at low temperature we performed a molecular 
dynamics simulation in which a superionic water sample in 
a 48-atom cell at a pressure of approximately 29 Mbar was 
quenched from a temperature of 5000 K to zero temperature 
over a period of 2 ps. In this simulation, as with all simulations 
described later, we used the VASP density functional theory 
code j2lTl . with pseudopotentials of the projector-augmented 
wave type J22ll . a cutoff for the expansion of the plane wave 
basis set for the wavefunctions of at least 1360 eV, and the 
PBE exchange-correlation functional B23I1 . The simulation 
used an N VT ensemble controlled by a Nose-Hoover thermo- 
stat in which the cell vectors were not allowed to change. In 
the resulting structure, the oxygen atoms retained their bcc lat- 
tice structure, while the hydrogen atoms all ended up close to 
the octahedral lattice sites, while none was close to a tetrahe- 
dral lattice site. This provides motivation for a more thorough 
investigation of this class of crystalline structures with octa- 
hedral hydrogen occupation. 

The bcc oxygen lattice provides four tetrahedral sites sur- 
rounding each oxygen and lying at the midpoint of the nearest- 
neighbor pairs, and six octahedral sites lying at the midpoint 
of second-nearest-neighbor pairs. This gives only one way 
to occupy the tetrahedral sites with two hydrogen atoms per 
oxygen, but many ways to occupy the octahedral sites with 
two hydrogen atoms per oxygen. On the assumption that the 
stablest configuration was likely to have reasonably high sym- 
metry, we systematically studied each possible configuration 
of hydrogen atoms for the three smallest possible unit cells: 
the six-atom lxlxl unit cell, the twelve-atom 2x1x1 
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FIG. 1: Cmcm ice structures where the large and small spheres de- 
note the oxygen and hydrogen atoms respectively. The orthorhombic 
unit cell with 12 atoms has been doubled in a direction and shifted 
so that an oxygen atom is at the origin. 



unit cell and the twelve-atom y/2 x y2 x 1 unit cell. This 
resulted, after the removal of equivalent structures, in three 
configurations for the lxlxl unit cell, six configurations 
for the 2x1x1 unit cell and also six configurations for the 
y/2 x v2 x 1 unit cell. Larger unit cells were found to have im- 
practically large numbers of inequivalent configurations. For 
each configuration, we performed a full geometry optimiza- 
tion at a pressure of 29 Mbar. 

The most stable configuration consists of a y2 x v2 x 1 
orthorhombic cell with 12 atoms and is shown in Fig. Q] Its 
reduced coordinates given in Tab. U The structure has Cmcm 
symmetry and may also be represented in 6-atom monoclinic 
unit cell with the vectors a m =(a + b)/2, b m =(a — b)/2, c m =c. 
This structure is one of a relatively small number of con- 
figurations which can be formed by filling two out of every 
three octahedral lattice sites with hydrogen atoms in such a 
way that no two vacancies are immediately adjacent. Dur- 
ing the structural relaxation, the hydrogen atoms that connect 
two oxygen atoms along the a direction move up along the 
b direction (y = 0.25 — > 0.34) towards the third vacant hy- 
drogen site. This elongates their bonds with the nearest oxy- 
gen atoms and introduces the kinks into O-H chains that are 
shown in Fig. Q] The structural relaxation of the next three 
lowest-enthalpy structures yielded enthalpies which are 0.51, 
0.59, and 1.1 eV per H2O unit higher. We also doubled and 
quadrupled the unit cell and re-relaxed structure starting with 
distorted hydrogen positions but no structure with lower en- 
thalpy was found. 

The orthorhombic structure with Cmcm symmetry is com- 
mon among high pressure materials including phase II of 
molecular hydrogen [24], AB and AB2 compounds, and the 
post-perovskite phase of ABO3. A Cmcm structure Il25ll with 
a four-atom unit cell was also proposed for the e phase of 
molecular oxygen but experiments later determined a differ- 
ent structure with C2/m symmetry lEill . 

Prompted by the referee, we also performed lattice dynam- 




FIG. 2: Structures of the ice X, Pbcm, and new Pbca and Cmcm 
phases are shown from top to bottom. The large and small spheres 
denote the O and H atoms respectively, while the thin lines denote the 
unit cells. The ice X to Pbcm transition is a displacement of atomic 
layers. In Pbca, the H atoms are squeezed out of midpoint between 
nearest O atoms. In Cmcm, the H atoms occupy mid-points between 
next-nearest O atoms. 



ics calculations with Abinit 12711 and VASP codes. This re- 
vealed a dynamic instability in the Pbcm structure. Fig. [3] 
illustrates that the mode (1/2,0,0) becomes unstable at 7.6 
Mbar. The structural relaxation in a 2 x 1 x 1 unit cell with 24 
atoms gave rise to a new intermediate structure of Pbca sym- 
metry (Tab. H} that precedes the transition to the Cmcm phase 
in pressure. The resulting sequence is shown in Fig. [2] In the 
Pbca structure, the hydrogen atoms are squeezed out of the 
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TABLE I: Structural parameters of different ice phases in orthorhom- 
bic unit cells. The last three lines specify the Wyckoff positions and 
the reduced coordinates of the atoms. The remaining positions fol- 
low from symmetry operations. 
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FIG. 4: Left: Electronic band structure for the Cmcm phase in a 
monoclinic unit cell at 22 Mbar. The arrows indicate where bands 
cross the Fermi level underlining the metallic nature of this phase. 
Right: Comparison of Pnma and Cmcm band structures at 12.5 Mbar 
in the orthorhombic unit cell illustrate the Peierls instability. 




0.2 0.3 
Phonon wavevector q,. 



FIG. 3: Dispersion curve of the Pbcm phonon band with the lowest 
energy along the a direction, exhibiting an instability at 7.6 Mbar. 



midpoint between two nearest oxygen atoms. However, they 
still reside near their tetrahedral sites. The distortion of the 
H positions occurs in alternating directions in the two hydro- 
gen bonded networks, which is accommodated by the unit cell 
doubling in a direction. 

The Cmcm structure is distinct from the other high-pressure 
ice phases in several ways. While ices VII, VIII and X as 
well as the Pbcm and Pbca structures all consist of two inter- 
penetrating hydrogen bonded networks, the Cmcm structure 
consists of corrugated sheets of H and O atoms (Figs. |T]&|2j. 
While in lower-pressure structures the hydrogen atoms occupy 
sites between nearest pairs of oxygen atoms, in the Cmcm 
structure they are squeezed into sites between second-nearest 
pairs of oxygen atoms. It should be noted that the characteri- 
zation of crystalline solids in terms of chemical bonds makes 
less sense at very high pressure because the structures are in- 
creasingly dominated by the mutual repulsion of atomic cores 
rather than the formation of favorable electronic bonds be- 
tween atoms. The repulsive forces that keep the atoms at their 



lattice sites increase strongly with pressure, however. 

Band structure calculations show that the Cmcm phase is 
metallic for all pressures under consideration. Fig. [4] shows 
that the band which was the conduction band in the lower- 
pressure structures, now dips below the Fermi level a the A 
point of the monoclinic cell and becomes partially occupied. 
Similarly a valence bands becomes unoccupied near V. 

However, careful structural relaxations in the pressure 
range of 7.6-15.5 Mbar reveal the existence of a Peierls in- 
stability that opens a band gap (Fig. [U by shifting hydrogen 
atom slightly away from the mid-point between near-nearest 
oxygen atoms. This distortion is small and fractional coordi- 
nates change by less than 2%. It lowers the enthalpy by just 
7 meV per molecule. The resulting structure has Pnma sym- 
metry with 12 atoms in an orthorhombic unit cell. Table U 
lists the structural parameters in the conventional coordinate 
setting for Pnma that differs from that for Cmcm. 

At 15.5 Mbar the band gap in the Pnma structure closes 
and the mechanism for the Peierls instability disappears. This 
pressure marks the transition from the Pnma to the Cmcm 
structure. Since all lower pressure ice phases are insulating, 
the transition to Cmcm at 15.5 Mbar marks the insulator-to- 
metal transition in water ice. 

We now compare the enthalpy, H = E + PV, of the ice 
X, Pbcm, Pbca, Pnma, and Cmcm structure after having op- 
timized the geometry in constant-pressure variable-cell simu- 
lations from 5 to 50 Mbar. Figure [5] shows the enthalpy dif- 
ference of the static lattice relative to the Pbcm structure. The 
ice X structure transforms into the Pbcm structure for pres- 
sures above approximately 3 Mbar, consistent with the results 
of Benoit et al. B12I1 . The Pbcm structure transforms into the 
Pbca phase at 7.6 Mbar consistent with our lattice dynamics 
calculations. 
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According to static lattice calculations, the Pbca phase 
would transform into Cmcm structure at 17.7 Mbar. To test the 
importance of zero point motion, we perform lattice dynamics 
calculation in the harmonic approximation with a 4 x 4 x 4 
q point grid. The Cmcm and Prima phases have significantly 
lower zero point energy than the Pbca phase (e.g. 75 meV per 
molecule at 12 Mbar) because the hydrogen atoms are less 
constrained at their octahedral sites between the more distant 
next-nearest oxygen atoms, leading to softer phonon modes. 
The resulting reduction in zero point energy means that the 
Pbca structure transforms into the insulting Pnma structure 
at 12.5 Mbar before this structure changes to Cmcm at 15.5 
Mbar. 

Prompted by the referee, we also tested the validity of the 
PAW VASP pseudopotentials by performing full-potential all- 
electron calculation with the Exciting code B28I1 . We focused 
this test on the Pbcm-to-Cmcm transition that occurs at 13.1 
Mbar (Pig. |5j within the PBE approximation. We recalculated 
this transition pressure within the local density approximation 
(LDA) and obtained 11.5 Mbar. Such a pressure difference is 
not unexpected because LDA typically underestimates transi- 
tion pressures 12911 . We then determined the energy from all- 
electron calculation for the five geometries in each phase that 
we obtained with LDA PAW calculations from 9 to 14.3 Mbar. 
We fit the resulting equations for state for both phases and 
recalculated the transition pressure. The all-electron method 
shifted the resulting transition pressure by less than 0. 1 Mbar, 
which confirms that the PAW pseudopotentials are sufficiently 
accurate for the purpose for this study. 
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FIG. 5: Enthalpy differences of the ice X, Pbca, Pnma, and 
Cmcm structures from the Pbcm structure as a function of pressure. 
PBE 1 23] was used except where LDA is indicated. Results are for 
static lattices without zero point motion. 

We have predicted that water ice will attain novel crystal 
structures at megabar pressures and determined the following 
sequence of structural transformations. At 7.6 Mbar the Pbcm 
phase will transform into a Pbca phase, which then changes 
into an insulating Pnma structure at 12.5 Mbar. At 15.5 Mbar 



a insulator-to-metal transition leads to a structure with Cmcm 
symmetry. This last transition is expected to greatly increase 
reflectivity in water ice, which will make it easier detect with 
spectroscopic techniques in dynamic high pressure experi- 
ments. Ramp compression fgl [l^l and pre-compressed i30ll 
shock waves appear as most promising techniques. 
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